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I.  INTRODUCTION 


For  more  than  fifty  years  it  has  been  recognized  that  our  understanding 
of  turbulent  flows  is  very  incomplete.  A  quotation  attributed  to  Sir  Horace 
Lamb  in  1932  might  still  be  appropriate. 

I  am  an  old  man  now,  and  when  I  die  and  go  to  Heaven  there 
are  two  matters  on  which  I  hope  for  enlightenment.  One  is 
quantum  electrodynamics  and  the  other  is  the  turbulent 
motion  of  fluids.  And  about  the  former  I  am  rather 
optimistic. 

According  to  Hinze,^ 

Turbulent  fluid  motion  is  an  irregular  condition  of  flow  in 
which  the  various  quantities  show  a  random  variation  with 
time  and  space  coordinates  so  that  statistically  distinct 
average  values  can  be  discerned. 


We  are  all  familiar  with  some  of  the  differences  between  laminar  and 
turbulent  flows.  Usually,  higher  values  of  friction  drag  and  pressure  drop 
are  associated  with  turbulent  flows.  The  diffusion  rate  of  a  scalar  quantity 
is  usually  greater  in  a  turbulent  flow  than  in  a  laminar  flow  (increased  "mix¬ 
ing")  and  turbulent  flows  are  usually  noisier.  A  turbulent  boundary  layer  can 
usually  negotiate  a  more  extensive  region  of  unfavorable  pressure  gradient 
prior  to  separation  than  can  a  laminar  boundary  layer. 

The  subject  of  turbulence  has  absorbed  the  energies  of  countless  research 
workers  over  a  period  of  more  than  two  decades.  It  still  continues  to  be  an 
area  of  research  where  lack  of  complete  understanding  prevails.  A  listing  of 
excellent  references  on  the  subject  of  turbulence  is  given  in  the  Bibliography 
at  the  end  of  this  report.  Turbulent  motion  is  a  time-dependent  phenomenon 
and  thus,  the  unsteady  Navier-Stokes  equations  are  considered  to  govern  the 
time  history  of  a  fluid  particle  in  a  turbulent  flow.  Numerical  procedures 
are  available  for  solving  such  equations;  however,  it  is  almost  impossible  to 
completely  analyze  a  turbulent  flow  in  this  way.  The  main  problem  is  that  the 
time  and  space  scales  of  the  turbulent  motion  are  extremely  small.  A  grid 
fine  enough  to  resolve  the  small  scale  motions  of  turbulence  would  therefore 
require  an  immense  and  impracticable  number  of  points.  The  number  of  grid 
points  required  and  the  small  size  of  the  time  steps  puts  the  practical  compu¬ 
tation  of  turbulent  flows  by  this  means  outside  the  realm  of  possibility  for 


7.  J.  0.  Hinze,  Tuvbu Lenae.  2nd  Edition,  MoGvoD-Hill,  New  York,  1976. 
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present  computers.  Some  researchers  are  optimistic  that  by  the  turn  of  the 
century  computer  technology  will  have  advanced  to  where  turbulent  flow  calcu¬ 
lations  can  be  made  from  first  principles. 

The  main  thrust  of  present  day  research  in  computational  fluid  mechanics 
and  heat  transfer  in  turbulent  flows  is  through  the  time  averaged  Navier- 
Stokes  equations  which  historically  have  been  referred  to  as  the  Reynolds 
equations  of  motion  in  many  circles.  Time  averaging  the  equations  of  motion 
gives  rise  to  new  terms  which  can  be  interpreted  as  new  "apparent"  stress 
gradients  and  heat  flux  quantities  associated  with  the  turbulent  motion. 
These  new  quantities  must  be  related  to  the  mean  flow  variables  through  turbu¬ 
lence  models  which  introduce  further  assumptions  and  approximations.  Thus, 
this  attack  on  the  turbulent  flow  problem  through  solving  the  Reynolds  equa¬ 
tions  of  motion  does  not  follow  entirely  from  first  principles  since  addition¬ 
al  assumptions  must  be  made  to  "close"  tho  system  of  equations.  This  closure 
is  achieved  via  turbulence  models.  A  turbulence  model  consists  of  a  set  of 
differential  equations  and/or  algebraic  equations  and  associated  constants, 
the  solutions  of  which,  in  conjunction  with  the  equations  of  mean  motion, 
closely  simulate  the  averaged  character  of  real  turbulent  flows. 

Various  turbulence  models  have  been  proposed.  Zero-equation  models  are 
useful  in  engineering  applications  but  their  applicability  is  limited  to  near 
equilibrium  flows. ^  One-equation  models  increase  the  computational  work  and 
do  not  bring  improvement  in  universality  and  predictive  capability  that  would 
justify  their  use.  Two-equation  models  are  more  universal  and  less  empirical. 
Higher  order  Reynolds  stress  models  are  sophisticated  and  have  gained  less 
popularity  than  two-equation  models.  In  the  present  study  numerical  computa¬ 
tions  are  made  using  Chien's^  k-e  t_wo-equatfon  turbulence  model  which  is 
similar  to  that  of  Jones  and  Launder.**  ®  Calculations  are  extended  up  to  the 
wall  and  the  exact  values  of  the  dependent  variables  at  the  wall  are  used  as 
boundary  conditions. 


2.  A,  J.  Wadaoak,  "Simple  Turbulence  Models  and  Their  Applications  to 
boundary  Layer  Separation,  "  NASA  CR-3283,  May  1980. 

3.  'uei-Yuan  Chien,  "Predictions  of  Channel  and  Boundary- Layer  Flows  with  a 
L aJ-Reynolds-Number  Turbulence  Model,  "  AIAA  Journal,  Vol.  20,  January 
"9‘^2,  pp.  33-38. 

4.  W.  P.  Janes  and  B.  E.  Launder,  "The  Prediction  of  Laminarizatian  with  a 
Two-Equation  Model  of  Turbulence,  "  Int.  Journal  of  Heat  and  Mass  Trans¬ 
fer,  Vol.  15,  1972. 

5.  W.  P.  Jones  and  B.  E.  Launder,  "The  Calculation  of  Low-Reynolds-Number 
Phenomena  with  a  Two-Equation  Model  of  Turbulence,  "  Int.  Journal  of  Heat 
and  Mass  Transfer,  Vol.  16,  1973. 

6.  B.  E.  Launder,  C.  H.  Pridden,  and  B.  I.  Sharma,  "The  Calculation  of 
Turbulent  Boundary  Layers  on  Spinning  and  Curved  Surfaces,  "  Journal  of 
Fluids  Engineering,  March  1977,  pp.  231-239. 


Most  of  the  well-known  turbulence  models  lead  to  parabolic  systems  of 
partial  differential  equations  when  coupled  with  the  conservation  equations. 
^  large  variety  of  numerical  methods  for  solving  parabolic  partial  differen¬ 
tial  equations  have  been  used  to  calculate  boundary- layer  flows  with  various 
levels  of  success.  The  Hartree-Womersley  method^  for  solving  parabolic 
equations  has  been  employed  to  solve  the  equations  for  laminar  and  turbulent 
boundary  layer  flows.  The  method  treats  derivatives  in  the  transverse  direc¬ 
tion  as  ordinary  derivatives  and  expresses  the  streamwise  derivative  as  a 
finite-difference.  This  reduces  the  system  of  partial  differential  equations 
to  a  sequence  of  ordinary  differential  equations  to  be  solved  in  succession  as 
the  integration  proceeds  downstream  from  one  station  to  the  next.  The  two- 
point  boundary- value  problem  for  the  ordinary  differential  equations  can  be 
solved  by  the  shooting  method.  The  inconvenience  associated  with  handling  the 
two-point  boundary  conditions  has  contributed  to  an  apparent  movement  away 
from  this  method  to  more  conventional  finite-difference  procedures. 

Finite-difference  schemes,  ranging  from  simple  conventional  ones  to  more 
sophisticated  variants,  have  been  used  extensively.  Pletcher®  used  the 
DuFor^-Frankel  explicit  scheme  to  calculate  incompressible  and  compressible 
turbulent  boundary  layers.  The  stability  constraint  associated  with  the 
explicit  scheme  does  restrict  the  allowed  step  size  for  numerical  integration. 
Implicit  schemes  do  not  have  that  constraint.  Diagonal  dominance  however, 
plays  an  important  part  in  these  schemes.^  Examples  of  the  implicit  methods 
can  be  found  in  the  work  of  Patankar  and  Spalding, Harris^^  and  Blottner^^ 
mong  others.  Implicit  schemes  of  Crank-Nicolson's  type  have  been  applied  to 


7 .  D.  R.  Havtvee  and  J.  R.  Womevsley,  "A  Method  for  the  Numevical  or  Mechan¬ 
ical  Solution  of  Certain  Types  of  Partial  Differential  Equations,  "  Proc. 
Ro;-al  Soc.  London,  A161,  p.  SIS,  1937. 

8.  R.  H.  Pletcher,  "On  a  Finite-Difference  Solution  for  the  Constant 
Property  Turbulent  Boundary  Layer,  "  AIAA  Journal,  Vol.  7,  February  1969, 
pp.  305- Sll. 

9.  R.  S.  Hirsh  and  D.  H.  Rudy,  "The  Role  of  Diagonal  Dominance  and  Cell 
Reynolds  Number  in  Implicit  Methods  for  Fluid  Mechanics  Problems,  " 

■1  oumal  of  Computational  Physics,  16,  1974,  pp.  304-310. 

10.  S.  V.  Patankar  and  D.  B.  Spalding,  Heat  and  Mass  Transfer  in  Boundary 
Layers,  Intertext  Books,  London,  1 9 70 . 

11.  ■'■.  E.  Harris,  "Numerical  Solution  of  the  Equations  for  Compressible 
Laminar,  Transitional,  and  Turbulent  Boundary  Layers  and  Comparisons  mth 
Experimental  Data,  "  NASA  TR-R  368,  1971  . 

12.  F.  G.  Blottner,  "Finite  Difference  Methods  of  Solution  of  the  Boundary 
Layer  Equations,  "  AIAA  Journal,  Vol.  8,  1970,  pp.  193-205. 

11 


e  =  e  -  2u  )2  =  0  (y). 


In  Equation  (18)  the  symbol  e  is  used  instead  of  e  since  this  notation  pre¬ 
vails  in  the  literature.  The  Jones-Launder  model  has  been  proved  to  be  a 
powerful  tool  for  the  calculation  of  boundary  layers,  free  shear  layers,  and 
some  recirculating  flows.  Recently,  new  forms  of  the  model  have  been  proposed 
in  order  to  improve  flow  predictions.  ^ The  model  proposed  by  Chien^  is 
more  well  behaved  mathematically  near  a  solid  wall  and  is  utilized  in  this 
study. 

Chien  proposed  a  few  changes  which  included  replacing  the  term 


in  the  k  Equation  (17)  by 


This  term  has  the  advantage  of  being  linear  in  k.  Both  the  terms  (-r- — )2  and 

k  ^ 

—  are  important  only  very  close  to  the  wall  and  become  very  small  away  from 

y2 

the  wall.  Thus  the  solution  away  from  the  wall  is  not  adversely  affected  by 
the  inclusion  of  this  correction  term.  The  term  in  Jones-Launder  model 

3kl/? 

(— 5— )2  is  nonlinear  and  suffers  from  gruJ  resolution  and  stability  problems, 
oy 


Another  of  Chien' s  proposals  was  to  make  e  of  the  order  0(y2)  at  the 
wall.  In  order  to  maintain  consistency  at  the  wall,  a  new  term  is  added  to 
the  e  equation  which  is  given  by 

-  2v  —  exp  (-  j  y"^) . 

y^ 


42.  H.  G.  Hoffman,  "Impvoved  Form  of  ‘he  Lao  Reynolde  Number  k-£  Tuvhulenee 
Model,"  Phye.  Fluids,  Vol.  18,  Ma-'ch  1975,  pp.  309-312. 

43.  C.  H.  G.  Lam  and  K.  Rremhovst,  "A  Modified  Form  of  the  k-e  Model  for 
Predicting  Wall  Turbulence,  "  J ournal  of  Fluids  Engineering,  Vol.  103, 
September  1981,  pp.  456-460. 
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:her  than  This  is  considered  advantageous  from  the  computational 

int  of  view  since  c  has  a  finite  value  at  the  wall  but  e  goes  to  zero  at  the 
11.  To  prove  this,  let  us  expand  the  fluctuating  velocities  close  to  the 
11  in  Taylor  series. 


u'  =  u^(t)y  +  U2(t)y2  +  - 


v'  =  v^(t)y  +  V2(t)y2  +  -  (20) 

w'  =  w^(t)y  +  W2(t)y2  +  - 


bstituting  these  into  the  definitions  of  k  and  e. 


k  =  ^  (u"2  +  v"2  +  w"2) 


j  [(u^  +  vf  +  w2)y2  +  2(u^U2  +  v^V2  +  WjW2)y3  +  ....]  (21) 


k  =  0(y2) 


Id 


£  =  v[(4iil)2  +  (1^)2  +  (~)2  ] 


-  V 


[(u2  +  v^  +  w2)  +  4(u^U2  ''l''2 


e  =  0(1) 


2  (11  —  )2 
^  3y  ^ 


(U?  +  +  wf) 


Id 
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follows:  (i)  adding  terms  accounting  for  the  viscous  diffusion  of  k  and  e; 

(ii)  replacing  constants  Cj,  03  and  C  with  functions  of  the  turbulence 
Reynolds  number  Rj  where  ^ 


» 


(iii)  adding  two  new  terms  in  Equations  (14),  and  (15);  hence,  the  proposed 
equations  are: 


"T 


=  C  f  p 


M  M 


Jsi 

e 


(16) 

I 


Dk 
^  Dt 


3y 


[(u  + 


Ti  *1  / \ 2 

-37  ^  *  "T  (37’  - 


pe  -  2p  (--yy-)2 


(17) 


Oe 

ITt 


a 

3y 


C(y  +  -/)  ly]  +  Cjf^  Vj  -  c^f^P  k 


3y' 


where 


fl  =  1.0, 

f2  =  1  -  0.3  exp  (-Rf), 

f^  =  exp  [-2. 5/(1  +  0.02  R-]-)  ]. 


I 


•t 


Inclusion  of  the  last  terms  in  Equations  (17)  and  (18)  makes  equation  (18)  an 
equation  for  the  quantity 


e  =  e 


ak^/2  ^ 

2» 


(19) 
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determined  by  the  shear  stress  at  the  wall  the  kinematic  viscosity  v,  the 

fluid  density  p  and  the  distance  y  from  the  wall  is  invoked.  One  can  then 
compute  the  flow  properties  at  a  point  beyond  the  viscosity-affected  zone  and 
locate  there  the  grid  point  nearest  to  the  wall.  In  this  approach  the  bound¬ 
ary  conditions  are  not  applied  at  the  wall  but  near  the  wall  and  it  is  refer¬ 
red  to  as  the  wall-function  method.  As  a  result  the  computer  time  required 
for  the  solution  of  the  governing  equations  is  considerably  reduced,  since  the 
necessity  to  resolve  the  steep  gradients  of  the  dependent  variable  in  the 

viscous  sublayer  is  removed.  In  practice  the  location  y'*'  of  the  first  grid- 
point  is  taken  in  the  region 


40  <  y"^  <  100 


and  the  values  of  the  dependent  variables  are  calculated  from  (see 
Bibliography) 


u^  =  -^  tn  (Ey^),  E  =  9.0,  <  =  0.41, 


k^  =  ^  =  3.33, 


+  _  1  _  2.44 

^  ”  +  ~  + 

<y  y 

where 


and  u^  =  /t^/p  . 


Although  many  flows  obey  these  near-wall  "laws"  with  sufficient  accuracy, 
deviations  are  observed  in  flows  characterized  by  separation,  strong  accelera¬ 
tion,  steep  temperature  gradients,  etc.  Computation  of  these  flows  cannot  be 
based  on  near-wall  region  "universal  laws."  Therefore,  turbulence  models 
incorporating  the  important  influence  of  viscosity  very  close  to  the  wall 
should  be  utilized  so  that  calculations  can  be  extended  up  to  the  wall  itself 
and  the  exact  values  of  the  dependent  variables  at  the  wall  can  be  used  as 
boundary  conditions. 


Jones  and  Launder  extended  the  k-e  model  to  include  the  influence  of  vis¬ 
cosity  very  close  to  the  wall  by  modifying  Equations  (9),  (14),  and  (15)  as 


The  k-£  model  used  in  the  present  study  is  based  upon  the  one  developed 
by  Jones  and  Launder^*’^  where  k  is  the  turbulent  kinetic  energy  and  e  is  the 
turbulent  dissipation  rate.  Using  the  Cartesian-tensor  notations,  k  and  e  are 
defined  in  terms  of  the  velocity  fluctuations  and  their  gradients  as 

k 


3uf  aui 


The  turbulent  eddy  viscosity  Uj  is, 


k2 

^'T  =  £  • 


The  variables  k  and  e  are  determined  from  the  sol u*^ ion  of  the  following 
transport  equations: 


/^T  3k \  .  / 3u \2  __ 


V3y^ 


where  Cj ,  C2,  0|^,  are  empirical  constants  having  the  values:  c^  =  1.55,  C2 

=  2.0,  Oj^  =  1.0,  =  1.3.  Equation  (14)  is  an  approximation  to  the  exact 

transport  equation  for  k  which  is  derived  from  the  incompressible  Navier- 
Stokes  equations.  The  e  Equation  (15)  is  formulated  in  analogy  to  the  k  equa¬ 
tion  since  it  can  not  be  derived. 


Direct  effects  of  the  molecular  viscosity  on  turbulence  structure  are 
neglected  in  most  turbulence  models.  Viscous  effects  are  indeed  negligible 
throughout  most  of  the  flow,  but  become  important  in  the  immediate  vicinity  of 
a  wall.  An  approach  not  used  in  this  study  is  to  avoid  the  complications  of 
the  viscosity-dependent  region  adjacent  to  a  wall.  The  assumption  that  the 
mean  velocity  and  the  statistical  correlations  in  this  region  are  completely 


The  eddy  viscosity  for  the  outer  region  is  given  by 


^“T^outer  *'cp  '^wake 

where  F^j^e  “  J-max  ^max  <>f  C„fc  the  smaller  of  the  two  values. 

The  quantities  y^^^^  and  are  determined  from  the  function  F(y)  =  y  |a)( 

[1  -  exp(-  y  /A^)]  where  is  the  maximum  value  of  F(y)  and  yi^gj^  is  the 

value  of  y  at  which  it  occurs.  The  function  F|^ig|5(y)  is  the  Klebanoff 
intermittency  factor  given  by 


^k,eh<^)  ■ 

The  quantity  is  the  difference  between  the  maximum  and  minimum  total 

velocity  in  the  profile, 

“dif  “  +  "^Imx  ■  *  ""'lin 


and  for  boundary  layers,  the  minimum  is  defined  as  zero. 

The  outer  formulation  can  be  used  in  wakes  as  well  as  in  attached  and 
separated  boundary  layers.  For  free-shear  flow  regions  or  wakes,  the  Van 

Driest  damping  term  [exp(-  y‘*'/A‘^)]  is  neglected.  It  is  necessary  to  specify 

the  following  constants;  A'*'  =  26,  C^-p  =  1.6,  Cj^igtj  =  0.3,  =  0.25,  <  =  0.4 

and  K  =  0.0168.  This  type  of  simple  model  is  generally  inadequate  for  complex 
flows. One  model  that  has  been  used  successfully  to  predict  many  flows 
is  the  two-equation  k-e  model. 


40.  J.  J.  Govski,  T.  R.  Govindan,  and  B.  Lakshminavayana,  "Computation  of 
Three-Dimeneianal  Tunbulent  Shear  Flous  in  Comers,  "  AIAA  Paper  No.  83- 
1733,  July  1983. 

41.  P.  Van  Guliak,  "Apptiaation  of  the  k-e  Turbulence  Model  to  a  Turbulent 
Boundary  Layer  Solution  for  Flao  about  a  Spinning  Yawed  Projectile  at 
Mach  3,  "  Master's  Thesis,  University  of  Delaware,  June  1983. 
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The  strength  of  a  model  lies  in  a  combination  of  accuracy  and  general¬ 
ity.  One  should  not  expect  two-equation  models  to  predict  flows  any  more 
accurately  than  simpler  models;  but  simpler  models  need  more  extensive  adjust¬ 
ments  for  each  different  flow  condition.  Reynolds  stress  formulations  are 
still  under  development.  An  excellent  review  of  the  status  of  turbulence 
modeling  for  computational  aerodynamics  has  been  made  by  Marvin^^  and  the 
performance  of  various  models  are  discussed.  The  two-equation  models  appear 
to  perform  better  for  separated  flows  especially  in  the  recovering  regions 
downstream.  Based  on  the  above  considerations,  a  two-equation  (k-e)  model  was 
formulated  and  utilized  in  the  present  study.  Additionally,  Baldwin-Lomax 
algebraic  model  2**  was  used  to  perform  the  same  computations  for  comparison 
purposes.  Both  of  these  models  are  described  below. 

C.  Baldwin-Lomax  Algebraic  Model 

The  algebraic  eddy  viscosity  model  used  in  this  study  is  that  developed 
by  Baldwin  and  Lomax. 2“*  It  is  a  two-layer  model  in  which  an  eddy  viscosity  is 
calculated  for  an  inner  and  an  outer  region. 


Pf  "  ^^T^ inner 

V  <  V 

^  -^crossover 

Pf  “  ^^T^ outer 

y  '>  y 

^  -^crossover 

where  y  is  the  normal  distance  from  the  wall  and  Ycrossover  smallest 

value  of  y  at  which  values  from  the  inner  and  outer  formulas  are  equal.  The 
inner  region  is  based  on  the  Prandtl-Van  Driest  formulation 


(  )  .  =  P  jl.^  I  U)  I 

'^T'lnner  ^  '  ' 


(11) 


where 


£  =  <  y[l  -  exp(-  y^/A"^)]  , 


p  u  y 


and  |w|  is  the  magnitude  of  vorticity  given  by 


lull  =  ri—  —12  +  fix  iWx2  +  fiw  .  iMu-|l/2 

I  '  “-^ay  '3z  3y^  ^3x  3z'  ■' 


39.  J.  G.  Marvin,  "Turbulence  Modeling  for  Computational  Aerodynamics,  "  AIAA 
Paper  No.  82-0164,  January  1982. 


where  is  taken  equal  to  0.09.  Some  of  the  currently  available  two-equation 

models  are  the  Jones-Launder,**  Ng-Spalding, Saffman-Wi Icox, ^2  Wilcox- 
Traci^s  and  Wi Icox-Rubesin^**  models.  The  Jones-Launder  (k-e)  and  Wilcox- 
Rubesin  (k  -  models  are  the  most  popular  of  these  models.  Several 
researchers  have  made  computations  using  these  turbulence  models  and  compared 
the  calculated  results  with  experimental  data.  However,  the  comparisons  have 
not  revealed  any  of  the  models  as  definitely  superior  over  the  others. 

B.  Reynolds  Stress  Models 

These  are  models  which  do  not  assume  that  the  turbulent  shearing  stress 
is  proportional  to  the  rate  of  mean  strain  i.e. , 


-  p  u'v'  * 


(10) 


Transport  equations  are  developed  for  all  the  components  of  the  Reynolds 
stress  tensor.  Such  modeling  requires  the  solution  of  three  or  more  partial 
differential  equations.  The  Reynolds  stress  model  proposed  by  Daly  and 
Harlow^^  is  a  good  example  which  requires  the  simultaneous  solution  of  five 
transport  equations.  To  date  these  models  have  been  used  largely  as  turbu¬ 
lence  research  tools. Recently,  these  models  were  applied  to  compute  the 
incompressible  flow  in  a  duct^'^  and  to  simulate  the  near-wake  flow  of  a  flat 
plate. 


32.  P.  G.  Saffman  and  D.  C.  Wilcox,  "Turbulence  Model  Predictions  for  Turbu¬ 
lent  Boundary  Lauere."  AIAA  Journal,  Vol.  12,  No.  4,  19?4, 

33.  D.  C.  Wilcox  and  R.  M.  Traci,  "A  Complete  Model  of  Turbulence,  "  AIAA 
Paper  No.  76-351,  1976. 

34.  D.  C.  Wilcox  and  M.  W.  Rubeein,  "Progress  in  Turbulence  Modeling  for 
Complex  Flou  Fields  Including  Effects  of  Compressibility,  "  NASA  TP-1517, 
1980. 

35.  B.  J.  Daly  and  F.  R.  Harlow,  "Transport  Equations  in  Turbulence,  "  Physics 
of  Fluids,  No.  13,  p.  2634,  1970. 

36.  B.  E.  Launder,  G.  J.  Reece,  and  W.  Rodi,  "Progress  in  the  Development  of 
Reynolds  Stress  Closure,"  J.  Fluid  Mechanics,  Vol.  68,  1975. 

37.  V.  Reitman,  M.  Israeli,  a  d  M.  Wolfshtein,  '^Numerical  Solution  of  the 
Reynolds  Stress  Equations  in  a  Developing  Duct  Flow,  "  AIAA  Paper  No.  83- 
1883,  July  1983. 


38.  A  Sugavanam,  "Near-Wake  Computations  with  Reynolds  Stress  Models,"  AIAA 
Paper  No.  83-1696,  July  1983. 


(7) 


d  u 

and  Mj  no  longer  becomes  zero  when  -^  =  0.  Other  one-equation  models  have 

been  suggested,  the  most  notable  one  being  the  one  used  by  Bradshaw,  et  al.^^ 
The  turbulence  energy  equation  is  used  but  the  form  of  modeling  of  the  turbu¬ 
lent  transport  terms  deviate  somewhat  from  the  one  described  above. 

3.  Two-Equation  Models:  In  one-equation  models,  the  length  scale  is 
evaluated  by  an  algebraic  expression  dependent  upon  only  the  local  flow  param¬ 
eters.  Researchers  in  turbulent  flow  have  long  felt  that  the  length  scale  in 
turbulence  models  should  also  depend  on  the  upstream  history  of  the  flow  and 
not  just  the  local  flow  conditions.  An  obvious  way  to  provide  such  a  depen¬ 
dence  of  i  on  the  flow  is  to  develop  a  partial  differential  equation  for  the 
transport  of  i.  This  then  is  the  main  motivation  behind  the  two-equation 
models. 


The  two-equation  models  involve  an  additional  partial  differential 
equation  which  in  effect  provides  the  turbulence  length  scale.  Researchers 
have  experienced  better  success  by  solving  a  transport  equation  for  a  variable 
related  to  the  turbulence  length  scale  rather  than  the  length  scale  itself.** 
The  other  transport  equation  used  is  for  turbulent  kinetic  energy  and  is  the 
same  one  used  in  one-equation  models. 

One  of  the  most  popular  two-equation  models  is  the  k-e  model  of  Jones 
and  Launder. *"5  Here  e  is  a  turbulent  dissipation  rate  and  is  assumed  to  be 
related  to  other  model  parameters: 


e 


z  • 


(8) 


The  turbulent  eddy  viscosity  is  related  to  e  as: 


(9) 


31,  P.  Bvadehau,  D,  H.  Fervie,  and  N,  P.  Atwell,  "Calaulaticn  of  Boundary 

Layer  Development  Using  the  Turbulent  Energy  Equation,  "  Journal  of  Fluid 
Mechanics,  No.  28,  p.  693,  1967. 
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£  =  <y  ; 


K  =  .41. 


(4) 


It  is  well  known  that  the  turbulence  must  be  damped  out  very  near  the 
wall  in  the  viscous  sublayer.  Van  Driest^®  proposed  an  exponential  damping 
factor  as  suggested  by  analogy  with  the  way  in  which  velocity  fluctuations  are 
observed  to  decay  near  an  oscillating  flat  plate  in  laminar  flow.  This  mixing 
length  proposal  gives, 

Ji  =  Ky  (1  -  e"'^  (5) 


For  flows  with  heat  transfer  a  simple  model  can  similarly  be  developed  for  the 
apparent  turbulent  conductivity.  These  simple  models  have  been  modified  and 
used  with  considerable  success  to  compute  a  relatively  wide  range  of  turbulent 
flows. 10*24>30 

2.  One-Equation  Models.  Although  the  simple  algebraic  model  works 
remarkably  well  with  a  range  of  flow  situations,  it  has  the  shortcoming  of 

predicting  py  as  zero  whenever  *■ —  — - -  — 


This  is  not  true  under  all  conditions. 


For  example. 


center 


deficiency  can  be  corrected;  but  the  applicability  of  the  algebraic  models  is 


limited  to  near-equilibrium  flows.  Most  of  the 
include  regions  which  are  far  from  equilibrium. 


flows  in  the  real  world 
All  these  factors  provide 
for  vr  and  require  the 


motivation  for  considering  other  interpretations  for  py  and  require  the 
application  of  more  advanced  turbulence  models. 

One-equation  models  are  models  which  require  the  solution  of  a  trans¬ 
port  partial  differential  equation  for  the  turbulent  velocity  scale  to  evalu¬ 
ate  the  Reynolds  stress.  The  length  scale,  n  is  still  specified  algebraical¬ 
ly.  The  turbulent  velocity  scale  Vj  is  written  as  the  square  root  of  the 

turbulent  kinetic  energy  k  defined  as: 


k  =  2  (u' 


and  the  transport  equation  developed  for  k  is  usually  used.  Thus,  Uy  can  be 
written  as: 


29.  E.  R.  Van  Dnieet,  "On  Tunhulant  Flam  Neav  a  Wall,"  Journal  of  the  Aero¬ 
nautical  Sciences,  Vol.  23,  No.  11,  November  1956. 

30.  W.  C.  Reynolde,  "Computation  of  Turbulent  Flams,  "  Ann.  Rev.  Fluid  Meah., 
Vol.  8,  1976,  pp.  183-208. 
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Reynolds-averaged  Navier-Stokes  equations  are  used  in  this  study.  The  infor¬ 
mation  lost  in  the  averaging  process  is  supplied,  in  an  approximate  fashion, 
by  a  turbulence  model. 

Various  turbulence  models  have  been  proposed  over  the  years  ranging  from 
simple  algebraic  eddy-vi scosity  formulations  to  sophisticated  Reynolds  stress- 
equations  models.  Most  models  can  be  categorized  as  either  a  turbulent  eddy 
viscosity  model  or  a  Reynolds  stress  model.  For  simplicity  the  models  are 
described  as  applied  to  boundary  layer  flows. 

A.  Turbulent  Eddy  Viscosity  Models 

This  class  of  model  is  based  on  a  concept  originally  advanced  by 
Boussinesq  in  1877.  The  assumption  is  that  the  turbulent  shearing  stress  can 
be  related  to  the  rate  of  mean  strain  through  an  apparent  turbulent  viscosity: 


p  u' v‘  =  u 


T  3y  • 


Models  of  this  type  can  be  quite  simple  or  complex  depending  on  how  pj  is 

related  to  other  flow  variables.  A  brief  description  of  some  turbulent  vis¬ 
cosity  concepts  is  given  below  in  increasing  order  of  complexity. 

1.  Zero-Equation  (Algebraic)  Models.  One  of  the  most  successful  simple 
model  was  suggested  by  Prandtl  in  1925: 


^T  -  ^ 


where  t  is  a  mixing  length,  a  characteristic  length  scale  of  turbulence.  An 
excellent  explanation  of  the  origin  of  the  model  is  given  by  Launder  and 
Spalding  (see  Bibliography)  who  presented  it  in  a  form  analogous  to  that  for 

the  molecular  viscosity  as  given  by  kinetic  theory  of  gases.  As  a  result  of 

this  analog,  t  4^  can  be  interpreted  as  a  characteristic  velocity  of  turbu- 
oy 

lenct',  Vj;  and  I  can  be  regarded  as  a  mean  free  path  for  collision  of  globules 

of  tijid.  Thus,  pj  can  be  thought  of  as: 

Pj  =  pVjt  .  (3) 

The  mixing  length  i  is  specified  as  an  algebraic  function  of  local  flow  param¬ 
eters  in  the  simple  models.  Prandtl  observed  that  the  mixing  length  is  pro¬ 
portional  to  the  transverse  distance  i.e.. 


of  computing  such  complex  flows  on  various  geometric  shapes  are  available.  A 
strong  need,  however,  exists  for  a  general  turbulence  model  that  can  be  used 
to  compute  such  complex  flow  fields  of  practical  interest.  The  two-equation 
turbulence  model  has  less  empiricism  and  wider  applicability  to  a  class  of 
complex  fluid  problems  than  the  algebraic  model.  The  extension  of  the  Navier- 
Stokes  algorithm  to  this  model  would  thus  be  an  important  advance.  The  objec¬ 
tive  of  the  present  research,  therefore,  is  to  incorporate  into  a  time 
dependent,  thin-layer  Navier-Stokes  code,  a  two-equation  turbulence  model 
which  uses  the  same  implicit  algorithm  and  generalized  geometry  formulation. 

Numerical  computations  are  made  of  two  transonic  flows  (i)  attached  flow 
over  an  axi symmetric  projectile  and  (ii)  separated  turbulent  flow  over  an  axi- 
symmetric  bump  model.  Computations  are  also  made  for  the  same  flow  situations 
using  the  zero-equation  eddy  viscosity  turbulence  model.  Both  turbulence 
models  are  assessed  by  comparing  calculated  values  of  wall  pressure  distri¬ 
bution  and  profiles  of  velocity,  turbulent  kinetic  energy  and  Reynolds  shear 
stress  with  experimental  measurements. 


II.  TURBULENCE  MODELS 

The  computation  of  an  entire  turbulent  flow  field  by  direct  numerical 
solution  of  the  time  dependent  conservation  equations  is,  at  present,  impossi¬ 
ble  due  to  the  extremely  fine  grid-spacing  required  to  resolve  the  smallest 
significant  eddies  of  the  flow  and  the  extremely  small  allowable  time-step. 
This  approach  to  turbulent  flow  computations  requires  computers  with  storage 
and  speed  capabilities  far  beyond  those  currently  available.  The  usual  point 
of  departure  in  practical  applications  is  an  averaged  version  of  the  conserva¬ 
tion  equations.  The  averaging  process  introduces  new  unknown  variables,  which 
must  be  modeled  in  terms  of  other  quantities.  Averaged  equations  may  be 
derived  through  time-  or  mass-weighted  averaging  at  flow  field  points  or  by 
averaging  the  conservation  equations  over  space.  The  latter  technique,  known 
as  "subgrid  modelling"  or  "large  eddy  simulation,"  is  prohibitively  expensive 
for  solving  practical  problems. 

Large  eddy  simulation  is  a  powerful  research  tool^^"^®  and  falls  between 
the  direct  simulation  of  turbulent  flows  and  Reynolds-averaged  Navier-Stokes 
calculations  both  in  cost  and  accuracy.  Because  of  the  cost  involved,  this 
technique  is  not  used  for  engineering  flow  predictions  at  present.  This  kind 
of  approach  will  ultimately  provide  more  understanding  and  will  eventually 
guide  the  development  of  models  that  include  more  physical  information. 


2ff .  A.  Leonard,  "Panel  Dieaueeion:  Large  Eddy  Simulation  Techniques,  "  AIAA 
Paper  No.  83-1878-CP,  July  1983. 

27.  K.  Dang,  "Evaluation  of  Simple  Subgrid-Saale  Models  for  the  Numerical 
Simulation  of  Homogeneous  Isotropic  and  Anisotropic  Turbulence,  "  AIAA 
Paper  No.  83-1692,  July  1983. 

28.  P.  Moin,  "Probing  Turbulence  via  Large  Eddy  Simulation,  "  AIAA  Paper  No. 
84-0174,  January  1984. 
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modified  by  the  explicit  implicit-characteristic  algorithm  of  Reference  19. 
Viegas  et  al^^  have  used  this  calculation  procedure  to  study  various  shock- 
wave  turbulent  boundary-layer  interaction  flows  bot*'  at  transonic  and  super¬ 
sonic  speeds. 

Recently,  codes  have  been  developed^*^"^^  which  solve  the  compressible  set 
of  Reynolds-averaged  thin-layer  Navier-Stokes  equations  for  high  Reynolds 
number  flows.  Parabolized  Navier-Stokes^o  computational  technique  is  used  for 
the  computation  of  supersonic  flows  whereas  the  unsteady  Navier-Stokes 
codes2i’22  can  be  used  for  both  transonic  and  supersonic  computations.  The 
thin-layer  Navier-Stokes  equations  are  cast  in  strong  conservation  law  form. 
The  equation  formulation  allows  for  arbitrary  body  geometries  and  is  solved 
using  an  implicit,  approximately  factored,  finite  difference  scheme  by  Beam 
and  Wanning. The  turbulence  model  used  is  an  algebraic  two  layer  eddy  vis¬ 
cosity  model  reported  by  Baldwin  and  Lomax. 2**  Such  simple  models  contain  a 
large  amount  of  empiricism  which  tends  to  make  these  models  inadequate  for 
complex  turbulent  flows. 

Real  world  problems  such  as  the  transonic  turbulent  flow  over  a  projec¬ 
tile  are  complex  due  to  the  presence  of  shock  waves.  The  flow  field  is  char¬ 
acterized  by  shock  wave-boundary  layer,  viscous-inviscid  interactions,  and  the 
large  separated  flow  region  behind  the  projectile  base.  It  is  advantageous  to 
use  the  thin-layer  Navier-Stokes  computational  technique  described  above  in 
that  it  considers  these  interactions  in  a  fully  coupled  manner.  As  the  capa¬ 
bility  for  computing  more  complex  flows  expands,  the  need  to  develop  a  more 
general  turbulence  model  also  expands.  Navier-Stokes  codes  that  are  capable 


19.  R.  V.  MaoCovmaakj  "An  Efficient  Numenioal  Method  for  Solving  the  Time- 
Dependent  Compressible  Navier-Stokes  Equations  at  Hi^  Reynolds  Number,  " 
Computing  in  Applied  Medhanias,  AMD  Vol.  18,  ASME,  1976. 

20.  L.  B.  Sdhiff  and  J.  L.  Steger,  "Numerical  Simulation  of  Steady  Supersonic 
Viscous  Flaw,"  AIAA  Paper  No.  79-0130,  January  1979. 

21.  J.  L.  Steger,  "Implicit  Finite  Difference  Simulation  of  Flaw  About  Arbi¬ 
trary  Geometries  with  Application  to  Airfoils,  "  AIAA  Journal,  Vol.  16, 

No.  4,  July  1978,  pp.  679-686. 

22.  T.  H.  Pulliam  and  J.  L.  Steger,  "On  Implicit  Finite-Difference  Simula¬ 
tions  of  Three-Dimensional  Flow,  "  AIAA  Journal,  Vol.  18,  No.  2,  February 
1980,  pp.  159-167. 

23.  R.  Beam  and  R.  F.  Warming,  "An  Implicit  Factored  Scheme  for  the  Compres¬ 
sible  Navier-Stokes  Equations,  "  AIAA  Paper  No.  77-645,  June  1977. 

24.  B.  S.  Baldwin  and  H.  Lomax,  "Thin  Layer  Approximation  and  Algebraic  Model 
for  Separated  Turbulent  Flows,"  AIAA  Paper  No.  78-257,  January  1978. 

25.  M.  Visbal  and  D.  Knight,  "Evaluation  of  the  Baldwin-Lomax  Turbulence 
Model  for  Two-Dimensional  Shock  Wave  Boundary  Layer  Interactions,  "  AIAA 
Paper  No.  83-1697,  July  1983. 


turbulent  boundary  layers  by  Harris^*  and  Cebeci  et  al^^  among  others. 
Patankar  and  Spalding^°  obtained  a  finite-difference  scheme  by  expressing  each 
term  in  the  governing  equations  as  an  integrated  average  over  a  small  control 
volume  defined  by  the  grid.  Their  general  calculation  procedure  was  applied 
to  a  wide  variety  of  flows  with  considerable  success.  The  Thomas  algorithm 
for  solving  the  tridiagonal  system  of  equations  is  usually  employed.  The 
Keller  box  method^**  has  been  adapted  to  turbulent  boundary  layer  calculations 
by  Keller  and  Cebeci. 

The  method  developed  by  Keller  for  parabolic  problems  is  second  order 
accurate  on  an  arbitrary  nonuniform  grid  network.  The  box-scheme,  being  an 
implicit  method,  is  stable  with  no  restrictions  on  the  grid  size  in  the 
streamwise  direction.  This  method  requires  the  solution  of  a block-tridiag- 
onal  system  of  equations.  Blottner^**  proposed  a  Crank-Nicolson  scheme  using 
a  variable  grid  and  claims  that  his  scheme  has  the  same  accuracy  as  the  Keller 
box  scheme  and  is  more  efficient  for  parabolic  equations.  The  calculation 
methods  described  thus  far  are  applicable  to  boundary  layer  flows. 

For  many  flow  situations  where  shock  wave  turbulent  boundary  layer  inter¬ 
actions  are  important,  boundary  layer  techniques  are  inadequate.  For  such 
complex  flows  the  differential  equations  used  to  describe  the  mean  floiw  are 
the  Reynold-averaged  Navier-Stokes  equations.  Excellent  reviews  of  the 
closure  concepts  for  these  equations  can  be  found  in  the  Bibliography.  The 
mean  flow  and  the  turbulence  field  equations  are  solved  simultaneously  in 
References  15-17.  The  numerical  procedure  used  is  the  basic  second-order, 
predictor-corrector,  finite-difference,  time-splitting  method  of  McCormack,^® 


12.  T.  Cebeci,  A.  M.  0.  Smith,  and  G.  Mosinskie,  "Calculation  of  Compneseihle 
Adiabatic  Turbulent  Boundary  Layers,  "  AIAA  Journal,  Vol.  8,  November 
:870,  pp.  1974-1982. 

14.  G.  Blottner,  "Variable  Grid  Scheme  Applied  to  Turbulent  Boundary 
ayere,  "  Computer  Methods  in  Applied  Mechanics  and  Engineering,  4,  1974, 

np.  179-194. 

15.  T.  -7.  Coakley  and  J.  R.  Viegas,  "Turbulence  Modeling  of  Shock  Separated 
P nundary-Layer  Flows,  "  Paper  presented  at  the  Symposium  on  Turbulent 
Ph'ar  Flows,  University  Park,  PA,  April  1977. 

16.  .  R.  Viegas  and  T.  J.  Coakley,  "Nwnerical  Investigation  of  Turbulence 
Uodels  for  Shock-Separated  Boundary-Layer  Flows,"  AIAA  Journal,  Vol.  16, 
'Jo.  4,  April  1978. 

17.  -K  R.  Viegas  and  C.  C.  Horstman,  "Comparison  of  Multiequation  Turbulence 
"''iels  for  Several  Shock  Boundary-Layer  Interaction  Flows,"  AIAA  Journal, 
'PoJ.  17,  August  1979,  pp.  811-820. 

18.  W.  MacCormack,  "Numerical  Solution  of  the  Interaction  of  a  Shock  IVave 
O’  th  a  Laminar  Boundary  Layer,  "  Lecture  Notes  in  Physics,  Vol.  8, 
'rringer-Verlag,  1971,  pp.  151-162. 
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It  is  added  to  balance  the  molecular  diffusion  term  in  a  manner  exactly 
analogousto  the  added  term  in  the  k  equation.  The  exponential  ensures  that 
this  term  decays  rapidly  and  its  effect  is  felt  only  close  to  the  wall. 

Additionally,  the  definition  of  has  been  modified  to  include  the  damp¬ 
ing  effect  due  to  the  presence  of  the  solid  wall  In  the  manner  of  Van  Driest's 
proposal . 


=  0.09  [1  -  exp  (-  0.01  y^)]  . 


(23) 


In  the  Jones-Launder  model  this  coefficient  is  defined  in  terms  of  the  turbu¬ 
lent  Reynolds  number.  With  these  modifications  the  k-e  model  now  takes  the 
following  form: 
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where  is  given  by  Equation  (23)  and 


Cl  =  1.44 

C2  =  1.92  [1  -  0.3  exp  (-Rf)]  ,  Ry  =  ^ 


A  more  general  form  of  the  k-e  equations  than  described  above  is  given  by 
Launder,  Sharma  and  Pridden®  as 


26 


where  k  is  the  turbulent  kinetic  energy,  e  is  the  turbulent  dissipation  rate 
and  Uj.  is  the  turbulent  viscosity  and  are  given  as. 


k  =  Y  +  v'2  +  w'2) 


The  empirical  coefficients  in  Equations  (26)-(28)  are  given  below: 

C  =  0.09  exp  [-3. 4/(1  +  0.02R.)2] 

U  u 

k2 

Rt  = 

Cl  =  1.44 
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C2  =  1.92  [1  -  0.3  exp  (- 


=  2.0 


a  =1.3. 
€ 


Chien's  modifications  are  included  and  the  model  used  in  this  study  is  as 
fol lows: 


ar  ^  3)rJ  ^  ‘'t  IT  (ir  ir^ 


-  pc  -  2u~ 
■'n 


3  3f  3U-  9U' 

" ir  it)  ^  iir  (ir  ^  'aX^ 


-C2  P 


r//2 


where  is  the  distance  normal  to  the  surface  and  the  empirical  coefficients 
given  iiy  Equation  (23)  and  (29). 

III.  NAVIER-STOKES  COMPUTATIONAL  TECHNIQUE 
A.  Governing  Equations 

The  three  dimensional  thin- layer  Navier-Stokes  equations  are  presented 
and  are  then  reduced  to  the  axi symmetric  formulation.  To  enhance  numerical 
accuracy  and  efficiency,  coordinate  mappings  of  the  governing  equations  are 
employed.  This  brings  the  body  surface  onto  a  coordinate  surface  (body-fitted 
system)  and  clusters  grid  points  in  flow  field  regions  where  dependent  vari¬ 
ables  are  expected  to  undergo  rapid  changes  (boundary  layer  for  example).  In 
the  transformed  plane,  uniform  discretization  formulas  and  well-ordered 
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interior  grid  point  solution  algorithms  can  be  used.  Related  work  using 
transformed  equations  in  flow  field  applications  can  be  found  1n  References 
21,  44-46. 

The  governing  equations  are  transformed  using  the  following  general 
coordinate  transformations: 


4  =  C(x,y,z,t) 
n  =  n(x,y,z,t) 
C  =  c(x,y,z,t) 

T  =  t 


The  equations  are  written  in  strong  conservation  law  form  and  the  transforma¬ 
tion  retains  this  form.**^  The  resulting  transformed  three  dimensional  thin- 
layer  Navier-Stokes  equations  can  be  written  in  nondimensional  form  as 


3^q  + 


3-E  +  3  F  +  3,G  =  Re"^3 J 
C  n  C  C 


(32) 
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44.  G.  S,  Deiwevt,  "Numeviaal  Simulation  of  High  Reynolds  Number  Transonic 
Flows,  "  AIAA  Journal,  Vol.  13,  No.  10,  October  1975,  pp.  1354-1359. 

45.  P.  Kutler,  S.  R.  Chakravarthy ,  and  C.  K.  Lombard,  "Supersonic  Flow  Over 
Ablated  Nosetips  Using  an  Unsteady  Implicit  Numerical  Procedure,  "  AIAA 
Paper  78-213,  1978. 

46.  R.  W.  MacCormack  and  A.  J.  Paullay,  "The  Influence  of  the  Computational 
Mesh  on  Accuracy  for  Initial  Value  Problem  with  Discontinuous  or  Non¬ 
unique  Solutions,  "  Computers  and  Fluids,  Vol.  2,  1974,  pp.  339-361. 


47.  H.  Viviand,  "Conservative  Form  of  Gas  Dynamic  Equations,  "  La  Recherche 
Aerospatile,  No.  1,  January- February  1974,  pp.  65-68. 


» 


~  pV 

"  PW 

puV  +  n  p 

X 

puW  +  c^p 

PVV  +  HyP 

G  = 

PVW  +  CyP 

pwV  +  n^p 

pwW  +  t^p 

(e+p)V  -  n,.p 

(e+p)W  -  c^p 

s  = 


0 

{(4+5^+;^)  [0.5u  (u^+v^+w^)^+  KPr“^(Y  -  l)‘^(a^)^ 
+  ( w/3)  ( c,^u+  i;^v+  ii^w)  ( c^y  ^+  c^v^+  } 


The  velocities 
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represent  the  contravarient  velocity  components. 

The  Cartesian  velocity  components  (u,v,w)  are  retained  as  the  dependent 
variables  and  are  non-dimensionalized  with  respect  to  a^  (the  free  stream 
speed  of  sound).  Pressure  is  defined  as 


P  =  (y  -  l)[e  -  .5p(u2  +  v^  +  w^)] 


(34) 


where  y  is  the  ratio  of  specific  heats,  p  is  the  density  referenced  to  p^  and 
e  is  the  total  energy  referenced  to  p  a^.  The  additional  parameters  are  (<) 
the  coefficient  of  thermal  conduct! vityT  (u)  the  dynamic  viscosity,  (Re)  the 
Reynolds  number,  (Pr)  the  Prandtl  number,  and  (X)  which  through  the  Stokes 
hypothesis  is  (-2/3) w. 
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The  metric  terms  of  Equation  (32)  are  defined  from 


and 


"x  ' 

J(V4- 

J (x  y ,  -  y  X,) 
n  c  ns 

^z " 

J(y.z  -  z.y  ) 

54.  =  -x5  -  y5  -  z5 

t  T  X  -^T  y  T  z 

''<V5  - 

n.  =  -xn  -  yn  -  zn 
t  T  X  T  y  T  z 

=  -x5  -yc  -  zi; 

t  T  X  T  y  T  z 

(35) 


n 


X  y,Zc  - 
n-'c  C 


X,y  Zr 


where  J  is  the  Jacobian  of  the  transformation.  For  the  computation  of  turbu¬ 
lent  flows,  u  and  <  comprise  of  their  molecular  and  turbulent  counterparts. 
The  turbulent  contribution  wj  and  are  supplied  through  an  eddy  viscosity 

hypothesis  described  in  Section  II. 

The  "thin-layer"  approximation^^’ is  used  here  and  is  valid  for  high 
Reynolds  number  flows.  In  high  Reynolds  number  flows  one  usually  has  only 
enough  grid  points  to  resolve  viscous  terms  in  a  thin  layer  near  the  body  sur¬ 
face.  Essentially,  all  the  viscous  terms  in  the  coordinate  directions  (here 
taken  as  K  and  n)  along  the  body  surface  are  neglected  while  terms  in  the  c  or 
the  near  normal  direction  to  the  body  are  retained.  This  approximation  is 
used  because,  due  to  computer  speed  and  storage  limitations,  fine  grid  spacing 
can  only  be  provided  in  one  coordinate  direction  (usually  taken  as  the  near 
normal  direction)  and  the  grid  spacing  available  in  the  other  two  directions 
is  usually  too  coarse  to  resolve  the  viscous  terms. 

The  three  dimensional  set  of  equations  are  then  reduced  to  obtain  the 
azimuthal-invariant  or  n-invariant  equations'*®  by  making  use  of  two  restric¬ 
tions:  (1)  all  body  geometries  are  of  an  axi symmetric  type;  (2)  the  state 


48.  C.  J.  Nietubiaz,  T.  H.  Pulliam,  and  J .  L.  Stegev,  "Humeviaal  Solution  of 
the  Azimuthal^ Invariant  Thin-Layer  Navier-Stokee  Equations,  "  U.S.  Army 
Ballistic  Research  Laboratory,  Aberdeen  Proving  Ground,  Maryland,  ARBRL- 
TR-02227,  March  1980.  (AD  A085716)  (Also  see  AIAA  Journal,  Vol.  18,  No. 
12,  December  1980.) 
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variables  and  the  contravariant  velocities  do  not  vary  in  the  azimuthal  direc¬ 
tion.  Here  n  is  used  for  the  azimuthal  coordinate.  A  sketch  of  a  typical 

A 

axi symmetric  body  and  the  coordinate  system  is  shown  in  Figure  1.  The 

term  of  Equation  (32)  is  thus  reduced  to  a  source  term  of  the  n-invariant 
equations.  The  resulting  thin-layer  n-invariant  or  axisymmetric  Navier- Stokes 
equations  are  then  written  as 


3_^q  +  3^E  +  3^G  +  H  = 


-1  ^ 
Re  3^S 


(36) 


where 


0 

0 

pV  {R^  (U  -  ^^)  +  R^  (W  -  c^)} 
-  pV  R(V  -  n^  )  -  P/R 
0 


and  q,  E,  G  and  S  are  as  defined  in  Equation  (32).  The  metric  terms  of  Equa 
tion  (35)  are  modified  as: 


Equation  (36)  contains  only  two  spatial  derivatives  but  does  retain  all  three 
momentum  equations  thus  allowing  a  degree  of  generality  over  the  standard  axi- 
symmetric  equations.  In  particular,  the  circumferential  velocity  is  not 
assumed  to  be  zero.  This  allows  computations  for  spinning  projectiles  or 
swirl  flow  to  be  accomplished.  The  n-invariant  equations  have  been  used  in  a 
number  of  flow  field  applications**®’^^  and  is  utilized  in  the  present  study. 

B .  Numerical  Method 

An  implicit  approximate  factorization  finite-difference  scheme  in  delta 
form  is  used  as  described  by  Beam  and  Warming.^®  An  implicit  method  was 
chosen  because  it  permits  a  time  step  much  greater  than  that  allowedby  explic¬ 
it  schemes.  For  problems  in  which  the  transient  solution  is  of  no  interest, 
this  offers  the  possible  advantage  of  being  able  to  reach  the  steady  state 
solution  faster  than  existing  explicit  schemes. 

The  Beam-Warming  implicit  algorithm  has  been  used  in  various  applica¬ 
tions.  2°’^**  ’‘*®’®®  The  algorithm  can  be  first  or  second  order  accurate  in  time 
and  second  or  fourth  order  accurate  in  space.  The  equations  are  factored 
(spatially  split)  which  reduces  the  solution  process  to  one-dimensional  prob¬ 
lems  ft  a  given  time  level.  Central  difference  operators  are  employed  and 
the  algorithm  produces  block  tridiagonal  systems  for  each  space  coordinate. 
The  main  computational  work  is  contained  in  the  solution  of  these  block  tri¬ 
diagonal  systems  of  equations. 
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The  resulting  finite  difference  equations,  written  in  delta  form,  are 
(I  +  h6  a"  -  A.J'^V  A  J)(I  +  h6  C"  -  A  J 

-  hRe'^6^J“^M"j)  X  -  q")  =  -At{6^E"  +  (38) 

-  Re‘^6^S")  -  AtH"  -  A^J'^[(V^A^)2  +  (7^A^)2  ]Jq". 


Here  h  =  At  because  only  first  order  accuracy  in  the  time  differencing  is 
needed  for  the  steady  state  flows  which  are  considered  here.  This  choice 
corresponds  to  the  Euler  implicit  time  differencing.  The  6's  represent 
central  difference  operators,  A  and  V  are  forward  and  backward  difference 

^  p  JNp 

operators  respectively.  The  Jacobian  matrices  A  =  — ^  ,  C  =  — x  along  with  the 

3q  3q 

coefficient  matrix  M  obtained  from  the  local  time  linearization  of  S  are 
described  in  detail  in  Reference  22.  Fourth  order  explicit  (A^)  and  implicit 

(Aj)  numerical  dissipation  terms  are  incorporated  into  the  differencing  scheme 

to  damp  high  frequency  growth  and  thus  to  control  the  nonlinear  instabilities. 
A  typical  range  for  the  smoothing  coefficients  is  A^  =  (1  to  5)  At  with  Aj  = 

3a^.  Details  of  the  algorithm  and  the  finite  difference  equations  as  they 

apply  to  turbulence  field  equations  will  be  discussed  in  the  next  Section. 

C.  Initial  and  Boundary  Conditions 

Free  stream  values  are  used  as  the  initial  conditions  in  the  entire  flow 

field  domain  of  interest.  Unknown  values  of  q  on  the  boundaries  are  updated 

explicitly  and  Aq  =  q"  -  q"  is  set  to  zero,  leading  to  a  first  order  error 
in  time  at  the  boundaries. 

The  updated  values  of  q  are  obtained  along  the  body  surface  by  linear 
extrapolation  of  p,  U  and  V  in  inviscid  flow.  In  viscous  flow  pis  extrapolat¬ 
ed  and  U  =  V  =  0.  In  either  case  W  =  0  and  values  of  u,  v  and  w  are  obtained 
from  the  following  relation. 


ru 


V 

w 


0  ""^y^z 

-u  -  c,- 

°  (^x^z  -  ^zS^  ° 

V  -  n^ 

0  n,e 

-  c. 

L  y  X  y  X  j 

L  t  J 

(39) 
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For  an  ax i symmetric  body  spinning  with  angular  velocity  id,  one  would  impose 
the  condition  V  =  Rw  on  the  body  surface.  Pressure  on  the  body  surface  is 
obtained  by  numerically  integrating  the  following  equation. 


XX  Z  Z  ^  X  z  c. 


+  u3^c^  +  v(JRx^)  +  w3^(JRx^)] 


-  pU(c^u^  +  c^w^)  +  pV[Jx^R2  (V  -  n^)] 


Here  is  the  normal  pressure  gradient  at  the  body  surface.  Equation  (40) 
results  by  combining  the  three  transformed  momentum  equations. 

The  axis  singularity  is  handled  as  in  Reference  22  where  flow  variables 
are  not  required  at  the  axis  due  to  the  fact  that  the  required  flux  vectors 
are  zero  along  the  singularity.  At  the  far  field  boundary  free  stream  values 
are  specified.  At  the  downstream  boundary,  first-order  extrapolation  is  used 
for  >  1  while  extrapolation  and  the  condition  that  pressure  is  fixed  at 
are  used  for  M  <  1. 


IV.  SOLUTION  OF  k-e  EQUATIONS 
A.  Turbulence  Field  Equations 

The  k-e  equations  used  in  this  study  can  be  written  from  Section  II  as. 


Dk  _  3  p,^t  .  V  3k  T  , 

p  Dt  "  tjt:  w:^  ^t  ax: 

J  k  ,1  J  J  1 


-  pe  -  2u  — 


3u.  3u. 


Dt  3X 


(  .  ^^0  ■3X.^  l*^t  k  3X.  ^SX.  3xJ 


.  e^  P  e  -y*/2 
-  C2  P  r - 2u  —  e 

•^n 
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^here  is  the  distance  normal  to  the  surface.  The  turbulent  kinetic  energy, 
:  and  the  turbulent  dissipation  rate,  e  are  given  as 


k  =  I  (u'2  +  v'2  +  w'2) 


(43) 


3U'. 


(44) 


and  the  turbulent  viscosity  is  related  to  k  and  e  by. 


r  k2 

“t  ■  ~ 


(45) 


The  empirical  coefficients  in  Equations  (41),  (42)  and  (45)  are  given  below: 


=  0.09  [1  -  exp  (-O.Oly"")] 


k2 

R.  =  — 
t  ve 


c^  =  1.44 


C2  =  1.92  [1  -  0.3  exp  (-  R^2)3 


(46) 


<^3 


=  2.0 


"k 


=  1.0 


a  =1.3. 
e 


Expanding  and  using  the  continuity  equation,  we  can  write  Equations  (41) 
and  (42)  in  conservation  form  as 


» 


5  =  5(x,y,z,t) 
n  =  n(x,y,2,t) 
;  =  c(x,y,z,t) 
T  =  t  . 


(52) 


uation  (51)  can  be  written  in  transformed  coordinates  while  still  retaining 
e  conservation  form**^  as 


► 


3q  3E  3F  3G  ^ 
3t  ^  ^  ^  ^ 


(53) 
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r  1  ^  X,  2  ^  2  2,  3e  T  ^t  e  r  /  2  2 
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The  velocities  1),  V  and  W  are  the  contravari ant  velocities  defined  in  Equation 
'33)  and  J  is  the  Jacobian  of  the  transformation.  The  thin-layer  approxima¬ 
tion  described  previously  is  used  here. 

Based  on  the  assumptions  for  axi symmetric  flow  described  in  the  previous 
thapter,  the  n  derivative  term  in  Equation  (53)  drops  out  and  the  axisymmetric 
set  of  turbulence  field  equations  can  be  written  in  non-dimensional  form  as 


13  +  li  ^  i 

3t  3C 


vhere 


1 

^  =  J 


=  J 


J  ' y  z' 


,"t  ^  ^  3k 


,  t  ^  V  9e 


Re  r ,  1  1  1.,  1  2  .  2n 

-- J—  [(C,  *  tj,  *  «  w^) 

/  \  vn  pc  n  -1  2u  k 

+  (<;  u  +  c  V  +  5  w  ) -  -j-  -  Re  r - 

'  X  4  y  i  zc,'  J  J 


Re"^  Cl  f  [(C^  +  +  C^)(u^  +  v^  +  w^) 

^  J  k  X  y  z' ^  c,  c  C 

+  (c,  u^  +  V,  +  c  w  )^  -  Co  T  ^  -  Re’^  —  —  e'-^ 


rictly  based  on  local  information  and  results  in  such  drastic  change  in 
havi or. 

The  turbulent  kinetic  energy  profiles  at  the  same  selected  stations  are 
own  in  Figure  30.  The  profiles  at  two  stations  upstream  of  the  shock  are 
own  in  Figures  30(a)  and  (b)  and  compare  well  with  experimental  measure- 
nts.  Figures  30(c)  through  (h)  show  the  profiles  in  the  separated  region, 
e  peaks  are  well  predicted  by  the  two-equation  k-e  model  although  the  loca- 
on  of  the  peaks  are  slightly  underpredicted.  Comparison  of  the  profiles 
ry  close  to  the  wall  indicate  poor  predictions.  In  the  redeveloping  flow 
gion  after  reattachment  the  computed  profiles  with  the  k-e  model  are  in 
cel  lent  agreement  with  the  experimental  data.  This  is  where  we  have  seen 
lod  agreement  in  the  mean  velocity  profiles  as  well. 

The  location  of,  and  variation  in,  the  maximum  turbulent  shear  stress  are 
lown  in  Figure  31  and  32  respectively.  As  shown  in  Figure  31  the  locations 
e  well  predicted  by  the  k-e  model  and  are  in  close  agreement  with  the  exper- 
lental  observations  except  near  x/c  =  1,  i.e.,  where  the  aft  end  of  the  bump 
;  affixed  to  the  cylinder.  This  is  in  the  separated  region  and  the  disagree- 
int  is  even  more  clear  in  the  variation  of  the  maximum  turbulent  shear  stress 
lown  in  Figure  32.  Additionally,  Figure  31  clearly  shows  the  peak  of  the 
•ofiles  shifting  away  from  the  wall  from  x/c  =  .5  to  1. 

The  location  of  the  maximum  turbulent  kinetic  energy  and  its  streamwise 
jriation  are  shown  in  Figures  33  and  34,  respectively.  As  seen  in  Figure  33 
ie  location  of  the  peaks  further  shift  from  the  wall  from  x/c  =  .5  to  1.0  and 
len  falls  off  in  the  same  way  observed  by  the  experimental  data.  Although 
ie  trend  is  the  same,  the  calculations  underpredict  the  location  of  the  peaks 
-igure  33)  and  overpredict  the  values  of  the  maximum  turbulent  kinetic 
lergies  (Figure  34).  It  is  also  clear  from  these  figures  that  the  k-c  model 
'edictions  are  in  good  agreement  with  the  experimental  data  in  the  redevelop- 
ig  regi on  (x/c  =  1 .25) . 


VI.  CONCLUDING  REMARKS 


.  Objecti ves 

The  objectives  of  the  reported  research  were: 

(1)  formulate  the  k-e  turbulence  model  in  general  spatial  coordinates 
id  incorporate  it  into  a  compressible,  axi symmetric,  thin-layer  Navier-Stokes 
ide , 

(2)  apply  the  resulting  solver  to  two  transonic  flows  for  which  experi- 
?ntal  data  are  available, 

(i)  an  ogi ve-cyl i nder-boattai 1  projectile  configuration  at  M  =  .94 

Id  .97 


(ii)  an  axisymmetric  bump  configuration  at  M^  =  .875  which  involves 
3cal  shock  induces  separation 
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ence  of  a  strong  shock  wave  on  the  model.  The  pocket  of  supersonic  flow 
on  (M  >  1)  can  also  be  seen.  Figure  23  shows  the  pressure  contours  in  the 
I  field  near  the  bump.  As  can  be  seen  the  flow  over  the  front  portion  of 
model  expands.  The  shock  wave  can  be  clearly  seen  to  exist  on  the  aft 
ion  of  the  circular-arc  bump.  These  qualitative  features  are  predicted  by 
I  the  algebraic  and  k-e  models. 

Figure  24  shows  the  velocity  vectors  over  the  aft  end  of  the  bump  obtain- 
yith  the  algebraic  model.  Flow  separates  and  the  reverse  flow  region  can 
seen  in  this  figure.  To  show  the  separation  bubble  more  clearly  stream 
;tion  contours  are  plotted  in  Figure  25.  The  two-equation  k-e  model  pre- 
;ion  is  shown  in  Figure  26.  This  model  predicts  a  smaller  separated 
ion.  The  two-equation  models  generally  predict  poorly  in  the  separated 
ion  and  do  well  in  the  redeveloping  flow  region  after  reattachment. This 
It  will  be  discussed  in  a  later  section. 

Figure  27  is  a  plot  of  the  surface  pressure  distribution  as  a  function  of 
longitudinal  position.  The  surface  pressure  is  referenced  to  the  total 
ssure.  The  longitudinal  position  in  this  plot,  and  plots  to  follow,  is 
erenced  to  the  leading  edge  of  the  bump  excluding  the  fairing  i.e.,  the 
ersecting  point  of  the  arc  of  the  bump  with  the  cylinder  (see  Figure  19). 
puted  results  are  obtained  with  both  turbulence  models  and  are  compared 
h  experiment.  The  position  of  the  shock  wave  is  well  predicted  by  both  the 
els;  however  there  is  a  small  disagreement  in  the  region  downstream  of  the 
ck.  The  largest  discrepancy  is  about  15%  and  could  partly  be  due  to  the 

ge  grid  spacings  used  in  the  redeveloping  flow  region. 

Development  of  the  mean  velocity,  turbulent  shear  stress  and  turbulent 
etic  energy  profiles  over  the  aft  portion  and  just  downstream  of  the  bump 
shown  in  Figures  28,  29  and  30  respectively.  The  mean  velocity  profiles 
shown  in  Figures  28(a)  -  28(f).  Flow  separation  occurs  as  shown  in  Figure 
a).  Figures  28(a)  through  28(c)  show  the  mean  velocity  profiles  in  the 
arated  region.  As  pointed  out  earlier,  the  k-e  model  predicts  a  thin 

ersed  flow  region.  It  is  especially  true  at  the  stations  selected  in  the 
arated  region.  Elsewhere  in  the  separated  region,  however,  it  is  not 

ignificant  (see  Figure  26).  Although  a  thicker  separated  region  ispredict- 
by  the  algebraic  model,  the  profiles  are  poorly  predicted  by  both  the 
els.  Away  from  the  wall,  k-e  model  calculations  show  better  agreement  with 
experimental  data.  Poor  predictions  for  both  turbulence  models  can  be 
erved  in  Figure  28(d)  at  the  station  just  upstream  of  reattachment.  The 
evelopment  of  the  flow  after  reattachment  is  shown  in  Figures  28(e)  and 
f).  Here  the  k-e  model  produced  a  solution  that  more  closely  represents 
experimental  data  than  did  the  algebraic  model. 

Figure  29  shows  the  turbulent  shear  stress  at  selected  streamwise  sta¬ 
rs.  As  evident  from  this  figure,  the  k-e  model  predicts  the  turbulent 
ar  stress  profiles  which  are  in  close  qualitative  agreement  with  the 
erimental  data.  The  peaks  are  not  as  well  predicted  however.  Additional- 
the  location  of  the  peaks  shifts  further  away  from  the  wall  just  as  deter- 
ed  experimentally  and  k-e  model  successfully  predicts  the  rate  of  peak 
placement  as  shown  in  Figures  29(a)  through  29(g).  The  algebraic  model  on 
other  hand  predicts  sharp  increase  or  decrease  in  the  turbulent  shear 
ess  as  seen  in  Figure  29(c)-(e)  and  29(f)-(h).  It  grossly  underpredicts 
turbulent  shear  stress  in  Figures  29(d)  and  29(g).  The  algebraic  model  is 
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Bachalo  and  Johnson.  The  data  was  obtained  in  the  NASA  Ames  2x2-Foot 

Transonic  Wind  Tunnel  using  a  LDV  technique.  Calculated  results  using  both 
the  algebraic  and  the  k-e  models  are  compared  with  these  experimental  data. 

A  schematic  diagram  of  the  model  and  its  associated  flow  field  is  shown 
in  Figure  19.  The  model  consists  of  an  annular  circular-arc  bump  affixed  to  a 
thin-walled  cylinder  of  outer  diameter  15.2  cm.  The  bump  has  a  thickness  of 
1.9  cm  and  a  chord  length  of  20.3  cm.  Its  leading  edge  is  joined  to  the 

cylinder  by  a  smooth  circular  arc  of  radius  18.3  cm  that  is  tangent  to  the 

cylinder  at  3.33  cm  upstream  and  to  the  bump  at  2.05  cm  downstream  of  the 
intersection  of  the  arc  of  the  bump  with  the  cylinder.  In  other  words,  a 
fairing  is  used  in  the  leading  edge  of  the  bump.  The  flow  field  contains  a 
separated  region  which  is  induced  by  a  shock  wave. 

The  computational  mesh  for  this  case  was  obtained  using  a  hyperbolic  grid 
generation  scheme.  The  grid  generated  this  way  is  orthogonal.  The  full 

grid  is  shown  in  Figure  20  whereas  Figure  21  shows  an  expanded  view  of  the 

grid  near  the  bump  model.  Most  of  the  grid  points  are  clustered  on  the  aft 

portion  and  just  downstream  of  the  circular-arc  bump  in  the  flow  direction. 
The  grid  points  in  the  normal  direction  were  exponentially  stretched  away  from 
the  wall.  The  first  point  was  taken  to  be  .000010  away  from  the  model  surface 

which  correspond  to  y"*"  of  about  0.5.  The  number  of  grid  points  used  was  78 
and  40  in  the  longitudinal  and  normal  directions,  respectively. 

The  upstream  boundary  conditions  were  prescribed  by  uniform  free  stream 
conditions.  First  order  extrapolation  was  used  at  the  downstream  boundary. 
The  no  slip  boundary  condition  was  used  on  the  wall  and  free  stream  conditions 
were  used  at  the  far  field  outer  boundary.  For  the  turbulence  variables  with 
the  k-e  two-equation  model,  k  and  e  were  set  to  zero  on  the  wall  and  at  the 
upstream  boundary.  Zero  derivatives  of  k  and  e  were  used  at  the  outer  and 
downstream  boundaries. 

Figures  22  and  23  show  the  qualitative  features  of  the  flow  field  near 
the  bump  model.  Figure  22  is  a  Mach  contour  plot  and  clearly  indicates  the 
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to  those  discussed  previously  for  M  =  .94.  The  peak  values  occur  at  y  ==  25. 
k  increases  over  the  boattail  corner  (X/D  =  5.05  to  5.36)  and  then  decreases 
over  the  boattail.  The  turbulent  dissipation  rate  profiles  are  shown  in 
Figure  13.  These  profiles  behave  better  than  k  profiles  in  the  region  outside 
the  edge  of  the  boundary  layer  and  drop  off  to  small  values  without  the 
presence  of  any  humps  in  the  profiles  in  that  region.  As  expected,  the  peaks 

in  c  profiles  occur  closer  to  the  wall  (y^  =  10)  than  those  of  the  k  profiles 

(y^  =  25). 

Figures  14  and  15  show  the  turbulent  eddy  viscosity  profiles  obtained 
with  the  algebraic  model  and  the  k-e  model,  respectively,  and  are  plotted  in 
physical  y  coordinate,  rises  to  its  peak  and  then  drops  off  sharply  over  a 

very  small  distance  from  the  surface.  The  magnitudes  of  at  each  of  these 

longitudinal  stations  differ  in  both  the  model  predictions  and  are  clearly 
shown  in  the  next  Figure  16.  Figure  16  is  plotted  in  the  law  of  the  wall 
coordinate  and  shows  the  variation  of  near  the  wall  more  clearly.  The 

profiles  with  k-e  model  have  sharper  peaks  compared  to  those  obtained  with  the 
algebraic  model.  Algebraic  model  predicts  sharp  increase  (X/D  =  5.61  to  6.19) 
and  decrease  (X/D  =  5.05  to  5.36)  in  whereas  k-e  model  predicts  rather 

gradual  change  since  it  takes  into  account  the  upstream  effects.  Comparison 
of  \i^  profiles  at  X/D  =  5.36  and  5.61  shows  poor  agreement  and  comparison  at 

the  other  three  stations  shows  good  agreement.  This  kind  of  a  disagreement  is 
local  and  may  not  have  a  large  overall  influence  on  the  results. 

Figure  17  shows  the  mean  velocity  profiles  at  the  same  longitudinal  sta¬ 
tions.  There  is  very  slight  difference  between  the  computed  results  obtained 
with  both  turbulence  models.  Comparison  of  the  calculated  profiles  have  been 
made  with  experimental  data  at  X/D  =  5.05,  5.36  and  5.61  and  the  comparison  in 
general  shows  good  agreement.  The  slight  difference  in  the  computed  results 

and  experimental  measurements  is  for  the  X/D  =  5.36  case.  This  profile  is 

only  .06  calibers  downstream  of  the  boattail  corner  and  is  in  the  vicinity  of 
severe  expansion.  The  experimental  data  was  reduced  using  wall  static  pres¬ 
sure  measurements.  The  greater  the  distance  from  the  wall, the  more  the  veloc¬ 
ity  data  may  be  in  error.  This  is  particularly  true  just  downstream  of  the 
expansion  corner  where  the  profile  may  extend  through  the  expansion  fan  with 
significantly  varying  static  pressures.  A  small  error  in  experimental 

measurements  thus  could  account  for  the  slight  difference.  The  computed  and 
experimental  surface  pressure  coefficient  are  again  shown  in  Figure  18  and 

compare  favorably. 


Separated  Flow  Over  an  Axi symmetric 


Numerical  computations  have  been  made  for  a  transonic  turbulent  separated 
flow  over  a  bump  model.  All  the  computed  results  shown  are  for  M  =  0.875,  a  = 
0°  and  Re  =  13.6  x  10^/m.  Experimental  measurements  of  the  mean  flow  quanti¬ 
ties  as  well  as  the  turbulence  variables  for  the  same  model  have  been  made  by 
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corner  and  X/D  =  5.36  is  just  after  the  boattail  corner.  Because  of  the 

severe  expansion  at  the  boattail  junction,  the  turbulence  kinetic  energy  is 
increased  by  a  factor  of  two  between  these  stations.  It  then  drops  off  on  the 
boattail  as  shown  by  the  profiles  at  stations  X/D  =  5.61  and  6.19.  The  humps 
in  these  profiles  are  believed  to  be  the  result  of  the  interaction  of  the 

shock  and  expansion  waves  with  the  turbulent  boundary  layer  and  occur  outside 

the  edge  of  the  boundary  layer.  The  peaks  in  the  k  profiles  occur  at  y^  =  20 
although  the  peak  is  moved  slightly  further  away  from  the  wall  near  the  boat- 
tail  corner  i.e.,  between  X/D  =  5.05  and  5.36.  As  shown  in  Figure  7  the 

turbulence  dissipation  rate  profiles  show  identical  behavior  for  the  same 
stations  with  the  exception  that  there  are  no  humps  present  in  the  region  out¬ 
side  the  edge  of  the  boundary  layer.  Additionally,  the  peaks  now  occur  closer 

to  the  wall  at  y^  =  10.  This  agrees  with  the  observed  behavior  of  the  peaks 
in  Reference  41. 

Turbulent  eddy  viscosities  are  found  from  k  and  e  with  the  two-equation 
model  and  algebraically  using  Baldwin  and  Lomax  model.  These  are  referenced 
to  the  molecular  viscosity  \i„  and  plotted  in  Figures  8  and  9  for  the  same 

longitudinal  positions  discussed  above.  Figure  8  shows  the  profiles 

obtained  with  the  algebraic  model.  The  profiles  have  rather  flat  peaks  and  go 
to  zero  outside  the  boundary  layer.  It  drops  off  sharply  in  magnitude  near 
the  boattail  corner  i.e.,  X/D  =  5.05  to  5.36  and  then  rises  sharply  on  the 
boattail  as  seen  by  the  profiles  between  X/D  =  5.61  and  6.19.  The  algebraic 
model  is  based  on  local  information  and  such  sharp  increase  or  decrease  in 

results.  The  profiles  obtained  with  the  k-e  model  on  the  other  hand  shows 

gradual  change  in  on  the  boattail  as  seen  in  Figure  9.  The  profiles  have 

sharper  peaks  and  then  fall  off  to  values  other  than  zero  outside  the  edge  of 
the  boundary  layer.  Although  k  and  e  profiles  drop  off  to  practically  zero, 
k^/e  does  not  drop  off  from  its  peak  value  monotonical ly  with  increasing  dis¬ 
tance  from  the  surface  and  results  in  non-zero  u^'s.  The  mean  flow  gradients 

outside  the  boundary  layer  are,  however,  exceedingly  small  and  these  u.’s  in 
no  way  adversely  affect  the  solution  of  the  mean  flow  quantities. 

Figure  10  shows  the  mean  velocity  profiles  at  the  same  selected  stations. 
Velocity  profiles  obtained  with  both  turbulence  models  compare  well  at  X/D  = 
3.42  and  6.19.  Experimental  data  is  available  at  the  other  three  stations  and 
are  used  for  comparison  with  the  calculations.  Both  models  predict  almost  the 
same  profile  at  X/D  =  5.05  and  comparison  with  experiment  is  good.  Just  down¬ 
stream  of  the  boattail  corner  i.e,,  at  X/D  =  5.36  and  5.61,  comparison  of  the 
k-E  calculations  with  experiment  are  in  better  agreement  than  the  algebraic 
model  predictions.  Figure  11  is  a  plot  of  the  surface  pressure  distribution 
as  a  function  of  the  longitudinal  position  over  the  projectile.  The  rapid 
expansion  at  the  ogive  and  boattail  junctions  is  apparent.  Computed  results 
obtained  with  both  models  are  compared  with  experiment  and  the  results  are  in 
good  agreement.  A  small  improvement  of  the  results  with  k-e  model  can  be  seen 
on  the  boattail. 

Results  are  now  presented  for  another  Mach  number,  M  =  .97  where  strong 
shock/boundary  layer  interactions  occur.  Figure  12  shows  the  turbulent 
kinetic  energy  profiles  at  various  longitudinal  positions.  These  look  similar 


Before  performing  the  computations,  the  flow  field  domain  of  interest 
must  be  discretized.  The  computational  grid  used  for  the  numerical  computa¬ 
tions  was  obtained  from  a  versatile  grid  generation  program  developed  by 
Steger,  et  al.^^  This  program  allows  arbitrary  grid  point  clustering  thus 
enabling  grid  points  to  be  clustered  near  the  body  surface  and  is  based  on  the 
elliptic  grid  generation  scheme  advocated  by  Thompson,  et  al.^°  In  this 
method  the  grid  in  the  physical  plane  is  defined  by  the  solution  of  a  Laplace 
or  a  Poisson  equation  and  the  generated  grid  is  not  orthogonal. 

The  full  grid  is  shown  in  Figure  4.  The  computational  domain  extended  to 
four  model  lengths  in  front,  four  model  lengths  in  the  normal  direction  and 
four  model  lengths  behind  the  projectile.  Such  an  extended  domain  is  used  to 
eliminate  the  possibility  of  any  wave  reflection  back  on  to  the  model.  The 
grid  consists  of  78  points  in  the  longitudinal  direction  and  40  points  in  the 
normal  direction.  An  expanded  view  of  the  grid  near  the  model  is  shown  in 
Figure  5.  The  dark  region  near  the  model  surface  results  from  clustering  of 
grid  points  which  are  needed  to  resolve  the  viscous  boundary  layer  region. 
The  grid  points  in  the  normal  direction  were  exponentially  stretched  away  from 
the  surface  with  a  minimum  spacing  at  the  wall  of  .00002  0.  This  spacing 
locates  at  least  two  to  three  points  within  the  laminar  sublayer.  Clustering 
in  the  longitudinal  direction  was  used  at  X/D  =  3.2  and  5.3,  the  ogive  and 
boattail  junction,  respectively,  where  appreciable  changes  in  the  flow  vari¬ 
ables  are  expected. 

The  projectile  base  was  modeled  as  an  extension  of  the  7.0°  boattail  for 
a  distance  of  two  calibers.  The  surface  line  was  then  turned  parallel  to  the 
model  axis  for  the  remainder  of  the  wake  region.  The  base  flow  is  thus 
modeled  as  an  extended  sting.  A  review  of  free-flight  shadowgraphs  for  pro¬ 
jectile  shapes  at  transonic  speeds  does  show  the  wake  flow  to  follow  near  the 
boattail  angle  for  a  distance  of  one  to  three  calibers  before  turning  parallel 
to  the  flow  direction. 

Results  are  first  presented  for  =  .94  and  a  =  0.  The  turbulence  quan¬ 
tities  k  and  e  obtained  with  the  two-equation  turbulence  model  are  shown  in 
Figures  6  and  7,  respectively,  for  selected  longitudinal  stations.  One  of  the 
longitudinal  stations  selected  is  near  the  ogive-cylinder  junction  and  the 
others  are  located  either  near  the  boattail  junction  or  on  the  boattail 
itself.  Note  that  the  station  X/D  =  6.19  is  on  the  extension  of  the  boattail. 
The  k-e  model  prediction  is  compared  with  that  of  the  algebraic  model  at  this 
station.  Figure  6  shows  the  turbulence  kinetic  energy  profiles  in  the  law  of 
the  wall  coordinate.  The  station  X/D  =  5.D5  is  in  front  of  the  boattail 
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At  the  upstream  boundary  k  and  e  are  set  to  zero  while  a  first-order  extrapo¬ 
lation  is  used  at  the  downstream  boundary. 


Coupling  with  Mean  Flow  Equations 


The  Equation  (38)  is  solved  first  by  the  method  described  in  Chapter  III 
for  the  mean  flow  quantities.  Next  the  turbulence  field  Equations  (63)  are 
solved  using  the  just  computed  mean  flow  quantities.  Solution  of  Equations 
(63)  give  k  and  e  and  Equation  (56)  is  then  used  to  compute  u^.  This  then 

becomes  the  input  in  the  solution  of  Equation  (38)  for  mean  flow  variables  and 
this  process  is  continued  at  each  time  step  until  steady  state  results  are 
achieved.  The  solution  procedure  of  the  turbulence  field  equations  lag  that 
of  the  mean  flow  equations  by  one  time  step. 


V.  RESULTS 

Numerical  computations  have  been  made  for  two  transonic  turbulent  flow 
cases:  (i)  attached  flow  over  an  axisymmetric  projectile;  and  (ii)  separated 
flow  over  an  axisymmetric  bump  model.  Both  the  algebraic  and  the  two-equation 
k-£  eddy  viscosity  turbulence  models  were  used.  Computed  results  are  present¬ 
ed  in  the  form  of  surface  pressure  plots,  velocity,  turbulent  kinetic  energy, 
turbulent  dissipation  rate  and  Reynolds  shear  stress  profiles.  Comparison 
with  experimental  data  has  been  made  to  assess  the  performance  of  both  turbu¬ 
lence  models. 

A.  Attached  Flow  over  an  Axisymmetric  Projectile 

The  transonic  flow  field  about  a  projectile  configuration  with  a  turbu¬ 
lent  boundary  layer  has  been  computed.  All  the  computed  results  shown  are  for 
a  =  0°,  Re  =  13  X  10®/m  and  M  =  0.94  and  0.97.  Numerical  results  are  compared 
with  experimental  measurements^^’ which  were  performed  for  the  same  shape  in 
the  NASA  Langley  Research  Center  8  foot  Transonic  Pressure  Tunnel. 

The  model  geometry  is  shown  in  Figure  3.  It  is  an  artillery  projectile 
consisting  of  a  secant-ogive  nose,  a  cylindrical  mid-section  and  a  7°  conical 
afterbody  or  boattail  of  half  a  caliber  (one  caliber  =  one  diameter). 
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Thus,  it  is  equivalent  to  the  central  differencing  and  has  an  additional 
smoothing  term.  A  is  a  smoothing  coefficient  inserted  into  this  numerical 
dissipation  term  and  can  be  varied  between  0  and  1.  A  =  1  would  correspond  to 
an  overall  first-order  accuracy  in  C.  A  typical  range  of  A  used  in  the  compu¬ 
tations  is  .01  to  .1.  The  smoothing  term  is  treated  explicitly. 

F.  Initial  and  Boundary  Conditions 

The  k  and  £  equations  are  marched  in  time  until  steady  state  results  are 
obtained  and  thus,  solved  as  an  initial-boundary  value  problem.  Initial  con¬ 
ditions  i.e.,  profiles  of  k  and  e  are  needed  initially  in  the  entire  flow 
field  region.  The  initial  conditions  can  be  arbitrary  but  it  may  take  longer 
time  to  get  the  converged  solution.  Therefore  more  realistic  profiles  of  k 
and  e  need  to  be  prescribed.  This  is  based  on  the  balance  of  each  of  the 
terms  in  the  equations.**^  An  example  of  the  balance  of  the  terms  in  the  k 
equation  is  reproduced  here  from  Reference  41, 

It  is  clear  from  Figure  2  that  large  gradients  in  the  turbulence  vari¬ 
ables  occur  very  near  the  wall  and  source  terms  are  dominant.  Convection 
terms  are  negligible  near  the  wall  and  the  largest  terms  are  the  production 
and  dissipation  terms.  Based  on  this  local  equilibrium,  we  equate 


production  =  dissipation 


to  obtain  the  initial  k  and  e  profiles.  The  turbulent  viscosity,  appears 

in  the  production  term  and  is  obtained  from  the  solution  with  an  algebraic 

eddy  viscosity  model.  The  above  assumption  works  well  for  attached  wall 
bounded  flows  and  is  poor  for  separated  or  free  shear  flows. 

Since  calculations  are  extended  up  to  the  wall,  it  is  easier  to  specify 

the  boundary  conditions  on  the  wall.  At  the  wall,  the  dependent  variables  are 

zero. 


k  =  e  =  0 


(70) 


In  the  far  field  which  lies  outside  the  edge  of  the  boundary  layer,  zero 
derivatives  of  k  and  e  are  used. 
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As  discussed  in  the  previous  chapter,  fourth  order  dissipation  terms  are 
usually  added  on  the  right  hand  side  of  Equation  (62)  to  help  control  the 
numerical  instability.  For  the  turbulence  field  variables  convective  terms 
often  dominate  the  diffusion  in  the  far  field  (away  from  the  wall)  and  can 
cause  convective  instability.  To  overcome  this  difficulty,  numerical 

9E 

smoothing  based  on  upwind  schemes^^’^^  is  used.  The  convection  term  for 

example  is  differenced  as 


(68) 


where  t  is  a  numerical  flux  given  by 


j+l/2  2  ^^j+1 


E.  +  E.,  lU.+U.J 
r  -  J  J  {on) 

E._i/2 - -J - 2  ^j-1' 


and  E  =  qU  . 


Substituting  these  numerical  fluxes,  the  right  hand  side  of  Equation  (68)  can 
be  simplified  and  Equation  (68)  can  be  written  as. 
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reduce  it  to  two  one-dimensional  operators.  With  the  source  term  retained, 
one  can  include  it  in  either  spatial  operators  or  in  both  or  factor  it  out 
altogether.  In  the  present  study,  we  have  included  it  in  the  c  operator  and 
the  factored  scheme  becomes, 


[I  +  At  6^a"][(I  -  At  d")  +  At(6^B"  -  6^c")]  Aq" 


-  At(6.E"  +  6  g"  -  ah")  +  At  S"  . 

s  4  4 


Expanding  the  factors  gives  Equation  (61)  back  plus  additional  higher  order 
terms  such  as 

A  ^  9  r  <1  ^  ^ 


At^A^A  6^B  Aq' 


For  steady  state  solutions  Aq  goes  to  zero  and  thus,  the  approximate  factor¬ 
ization  error  vanishes.  The  factored  form  Equation  (62)  has  reduced  the  two- 
dimensional  matrix  inversion  problem  to  two  one-dimensional  problems  which  can 
be  efficiently  solved. 

E.  Solution  Algorithm 

A  convenient  solution  algorithm  is  developed  for  Equation  (62)  with  the 
following  sequence. 

[(I  -  Ato")  +  At(A^B"  -  A^c")]  Aq^  =  RHS  (62)  (63a) 


[I  +  AtA^A®]  Aq”  =  AcT 


q  =  q  +  Aq 


(63b) 


(63c) 


where  RHS  (62)  is  the  right  hand  side  of  Equation  (62).  First,  equation  (63a) 

is  solved  for  A^  since  the  right  hand  side  is  known  at  the  old  time  step. 

A^  becomes  the  right  hand  side  of  Equation  (63b)  which  is  then  solved  for 

Aq”.  This  is  then  added  to  q^  to  give  q'^^  at  the  next  time  step  as  shown  in 
Equation  (63c). 
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Substituting  Equation  (59)  into  (58)  and  rearranging,  one  obtains 


[(I  -  At  d")  +  At(6^A"  +  -  6^c")]  (q'^^  -  q") 


-  At(5^E"  +  -  a^H")  +  At  S" 


(61) 


where  I  is  an  identity  matrix  and  6^  ,  6^  are  the  spatial  difference  '■*  'a- 
tors.  This  is  in  the  "delta"  form  since  we  are  solving  for  Aq"  =  q^^  -  q*^. 

D.  Approximate  Factorization 

Direct  inversion  of  the  block  matrix  on  the  left  hand  side  of  Equation 
(61)  is  a  formidable  task.  To  avoid  this  problem,  approximate  factorization 
of  the  left  hand  side  operator  of  Equation  (61)  is  frequently  used.  In  the 
absence  of  the  source  term,  one  can  approximately  factor  the  left  hand  side  to 


The  source  terms  S  associated  with  the  turbulence  field  variables  can  be  very 
large.  As  happens  near  the  wall,  the  source  terms  (production,  dissipation 
and  decay)  become  dominant  over  the  convection  and  diffusion  terms.  This  can 
result  in  a  very  stiff  algorithm  if  the  source  terms  are  treated  entirely  in 
an  explicit  manner. S'*  Thus,  they  are  treated  implicitly  as  shown  in  Equation 
(58). 


Equation  (58)  is  nonlinear  since  E,  G,  H  and  S  are  functions  of  the 
dependent  variable  q.  The  nonlinearity  can  be  removed  by  a  linearization 

procedure.  A  local  Taylor  expansion  about  q^  yields. 


.  E"  +  (||)"  (q"^l  -  q")  +  0(At2) 


gn+1  ^  gn  ^  ^^n+l  _  ^n^  ^ 


=  H"  +  (|M)n  (qH+l  _  qO)  ^  q(^^2) 


gn+l  ^  gn  ^  (||jn  ^^n+l  .  ^n^  ^  ^^^^2)^ 


Let  A  =  -  —  B  =  - --  C  =  —  D  -  — 

^  3q  ’  “  3q  ’  ^  9q  ’  ^  ■  9q 


These  Jacobian  matrices  are: 
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The  variables  in  Equation  (54)  are  made  dimensionless  as  follows: 


P  = 


p_ 
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00 
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(55) 


For  simplicity  the  'bars'  have  been  removed  from  the  variables  in  Equation 
(54).  The  turbulent  viscosity  in  nondimensional  form  becomes 
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(56) 


where 


Re 


p„a  D 
00  00 

00 


C.  Numerical  Method 


The  numerical  scheme  used  is  the  Beam-Warming^^  Euler  implicit  scheme. 
The  time  differencing  is 


q"^^  =  q"  +  At  (||)"^^  +  0(At2) 


(57) 


If  Equation  (54)  is  inserted  in  (57),  one  has 
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(3)  compare  the  predicted  mean  flow  and  turbulence  quantities  with 
experimental  data,  and 

(4)  compare  the  k-e  turbulence  model  with  an  algebraic  mixing  length 
turbulence  model  as  applied  to  above  problems. 

B.  Summary  of  Results 


The  thin  layer  form  of  the  compressible  Navi er-Stokes  equations  was 
solved  using  a  time  dependent,  implicit,  approximately  factored,  finite  dif¬ 
ference  scheme.  The  equations  were  marched  in  time  until  the  desired  steady 
state  results  were  achieved. 

(1)  For  the  computation  of  turbulent  flows,  the  turbulence  closure 
was  provided  with  Baldwi n-Lomax  algebraic  and  Chien's  k-e  two-equation  eddy 
viscosity  models.  The  k  and  e  equations  were  developed  in  the  general  spatial 
coordinates  and  incorporated  into  a  thin  layer,  time  dependent  Navier-Stokes 
code.  The  same  implicit  algorithm  that  simultaneously  solves  the  mean  flow 
equations  was  extended  to  solve  the  turbulence  field  equations  using  block 
tridiagonal  matrix  inversions.  Calculations  with  the  k-e  model  have  been 
extended  up  to  the  wall  and  the  exact  values  of  the  dependent  variables  at  the 
wall  have  been  used  as  boundary  conditions.  Very  small  grid  spacing  was 
utilized  close  to  the  wall  in  order  to  resolve  the  steep  gradients  of  the 
dependent  variables  observed  in  the  viscous  sublayer.  The  distance  of  the 

first  grid  point  from  the  wall  should  be  within  y^  c  1.25. 

(2)  The  transonic  flow  field  about  an  artillery  projectile  (secant- 
ogive,  cylinder,  boattail)  with  a  turbulent  boundary  layer  has  been  computed. 
The  computed  results  were  obtained  for  a  =  0°,  Re  =  13  x  10^/m  and  =  .94 
and  .97.  These  results  were  compared  with  the  available  experimental  measure¬ 
ments  of  the  mean  flow  quantities. 

Computed  results  show  the  turbulent  kinetic  energy,  dissipation  rate  and 
turbulent  eddy  viscosity  profiles.  The  velocity  profiles  and  the  surface 
pressure  distribution  have  been  obtained  with  both  the  algebraic  and  the  k-e 
turbulence  models  and  are  compared  to  experiment.  The  results  are  in  good 
agreement.  The  rapid  expansions  at  the  ogive  and  the  boattail  junctions  are 
well  predicted  by  both  models.  A  small  improvement  with  the  k-e  model  predic¬ 
tion  is  found  at  M  =  .94. 

00 

(3)  Numerical  computations  have  also  been  made  for  a  transonic  tur¬ 

bulent  flow  over  an  axisymmetric  bump  model  which  involves  shock  induced  sepa¬ 
ration.  The  computed  results  were  obtained  for  =  .875,  a  =  0°  and 

Re  =  13.6  X  10^/m.  The  computed  results  are  compared  with  the  experimental 
measurements  for  both  the  mean  flow  and  the  turbulence  quantities  which  are 
available  for  this  case. 

The  surface  pressure  distribution  and  the  contour  plots  of  Mach  number 
and  pressure  indicate  the  presence  of  a  strong  shock  wave.  The  position  of 
the  shock  wave  is  well  predicted  by  both  the  algebraic  and  the  k-e  turbulence 
models  and  compares  well  with  experiment.  Results  are  presented  showing  the 
development  of  the  mean  velocity,  turbulent  shear  stress  and  turbulent  kinetic 
energy  profiles  over  the  aft  portion  and  just  downstream  of  the  bump.  The 


54 


.F.  V.  IPJ  '."JL".'  -■ 


results  are  generally  in  good  agreement  with  the  experimental  data.  Predic¬ 
tions  by  both  turbulence  models  are  poor  in  the  separated  flow  region.  In  the 
redevelopment  region  downstream,  however,  k-e  model  prediction  is  in  better 
agreement  with  the  data.  The  k-e  model  successfully  predicts  the  location  and 
the  trend  in  the  peaks  of  the  turbulent  shear  stress  and  turbulent  kinetic 
energy  profiles.  The  algebraic  model  predicts  sharp  increase  and  decrease  in 
the  turbulent  shear  stress  which  is  physically  unrealistic. 

(4)  The  algebraic  model  is  based  on  local  information  and  predicts 
undesirable  sharp  increase  and  decrease  in  the  turbulent  shear  stress.  As 
expected,  the  k-e  model  avoids  this  since  length  scales  are  obtained  by  solv¬ 
ing  a  transport  equation.  Poor  comparison  between  the  predictions  by  both 
models  and  the  experiment  was  found  in  the  recirculating  region.  Some 

improvements  were  found  in  the  developing  regions  downstream  with  the  k-e 
model.  Where  the  mean  velocities  are  relatively  in  good  agreement  with  the 
experimental  results,  so  are  the  turbulent  shear  stress  and  kinetic  energies. 

C.  Recommendations 


The  major  difference  between  the  calculations  and  the  experiment  is  in 
the  separated  region.  It  is  in  this  region  that  the  turbulent  shear  stress 
and  kinetic  energy  are  under  predicted  by  the  k-e  model.  At  this  stage  it  is 
exceedingly  difficult  to  sort  out  the  discrepancies  between  computation  and 
experiment  that  arise  separately  from  turbulence  modeling  and  computational 
procedures.  This  situation  will  continue  until  grid  independent  computations 
can  be  achieved  and  numerical  smoothing  procedures  are  fully  understood. 

With  this  in  mind  one  can  only  speculate  for  improving  the  model  predic¬ 
tions.  A  look  at  the  balance  of  the  terms  in  the  k-equation  suggests  that  the 
balance  of  the  production  and  the  dissipation  must  occur  farther  away  from  the 
wall  in  order  to  produce  experimental ly  observed  peaks  in  turbulent  shear 
stress  and  kinetic  energy.  This  in  turn  implies  tuning  the  e  equation. 
Further  computational  investigation  is  recommended  in  this  regard.  The 
protuberance  configuration  is  a  good  case  for  testing  turbulence  models.  It 
is  recommended  that  computations  be  made  for  this  configuration  to  further 
validate  the  k-e  model.  Additionally,  the  wake  flow  or  the  base  flow  behind 
the  base  of  the  projectile  is  one  where  k-e  model  application  is  more  appro¬ 
priate  since  it  predicts  gradual  change  in  length  scales.  It  is  believed  that 
the  body  of  information  obtained  in  the  present  research  forms  a  sound  basis 
for  attacking  such  complex  flow  problems. 
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Figure  1.  Axisymmetric  Body  and  the  Coordinate  System 
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Figure  2.  Balance  of  Terms  in  the  k-equation  (Reference  41) 
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Figure  6.  Turbulent  Kinetic  Energy  Profiles,  M  =  .94,  a  =  0 
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Figure  9.  Turbulent  Viscosity  Profiles,  =  .94,  a  =  0  (k-E  Model 
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Figure  12.  Turbulent  Kinetic  Energy  Profiles,  M, 
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Figure  15.  Turbulent  Viscosity  Profiles,  =  .97,  a  =  0  (k-e  Model) 


Figure  16a.  Turbulent  Viscosity  Profiles,  M, 


0,  X/D  =  3.42 


y  (cm) 


0  5  10  15  20 

-uv/u»  «  10 

Figure  29d.  Turbulent  Shear  Stress  Profiles,  =  .875,  a  =  0,  x/c  =  .87! 


Figure  29e.  Turbulent  Shear  Stress  Profiles,  M 
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Figure  29b.  Turbulent  Shear  Stress  Profiles,  =  .875,  a  =  0,  x/c  =  .625 


Figure  29c.  Turbulent  Shear  Stress  Profiles,  M, 
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Figure  28g.  Mean  Velocity  Profiles,  M 
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Figure  28e.  Mean  Velocity  Profiles,  M„  =  . 
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Figure  28c.  Mean  Velocity  Profiles,  M, 
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Figure  28a.  Mean  Velocity  Profiles,  M, 
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Figure  25.  Stream  Function  Contours,  =  .875,  a  =  0  (Algebraic  Model) 


Figure  26.  Stream  Function  Contours,  M^  =  .875,  a  =  0  (k-e  Model) 
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Figure  24.  Velocity  Vectors, 
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Figure  19.  Schematic  Illustration  of  the  Bump  Model 
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Figure  29h.  Turbulent  Shear  Stress  Profiles,  =  .875,  a  =  0, 


k-e  MODEL 
ALGEBRAIC  MODEL 


O  EXPERIMENT 


E 

1.0 

>. 


■UvVuo**  *  10* 


Figure  29i.  Turbulent  Shear  Stress  Profiles,  M  =  .875,  a  =  0 
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Figure  30c.  Turbulent  Kinetic  Energy  Profiles, 
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Figure  30d.  Turbulent  Kinetic  Energy  Profiles,  =  .875,  o  = 


Figure  30e.  Turbulent  Kinetic  Energy  Profiles,  M 
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Figure  30h.  Turbulent  Kinetic  Energy  Profiles,  =  ,875, 
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Figure  30i.  Turbulent  Kinetic  Energy  Profiles,  M_  =  .875 
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Figure  30j.  Turbulent  Kinetic  Energy  Profiles,  =  .875,  a 
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Figure  31.  Location  of  Maximum  Turbulent  Shear  Stress,  M. 
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Figure  34.  Variation  of  Maximum  Turbul 
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